NW GOM ROV analyses

version: March 27, 2020

About this document

---

All analyses preformed with R version 3.6.1.

This is the code that accompanies the publication XXX. Here you will find all the code to repeat the statistical analyses performed for this manuscript. All of the accompanying data can be found on GitHub.

For this project we conducted 296 ROV transects (~0.1 km length) in mesophotic coral ecosystems (MCEs) across 5 banks in the nortwestern Gulf of Mexico (NW GOM). Downward-facting still photographs including scaled parallel laser points were captured along each transect. These photographs were analyzed for benthic percent cover using 50 randomly placed point per photograph.

If no biota were present under a particular point, the benthos was characterized as soft or hard bottom. Biota were identified to the lowest taxonomic level possible, given the image quality. Scleractinian corals were identified to species, cnidarians in the Class Alcyonacea and Order Antipatharia were identified to family level, sponges were identified to Class, and algae to Phylum. All other categories were categorized into bacterial mat, Chordata, echinoderm, hydroid, Mollusca, and Tunicata groups. Dead coral was marked as different from hard bottom substrate, and a miscellaneous category was used for unidentifiable biota. Points that could not be identified at all, such as those in shadows or sediment plumes, were labeled as “tape wand shadow” and were removed from statistical analyses.

Coral density was also calculated as individuals m^-2 across the entire transect image aresa using the scaled lasers to calculate image sizes.

If you download my entire accompanying github directory you should be able to re-run these analyses by following along with the code chunks in R Studio. If you download the code separtely or you are using this pipeline on your own data, you may need to change the working directory to where the associated files are housed (ie. setwd("~/path/to/directory/with/data")).

Basic setup of R environment

---

Loading required packages

For the following analyses we will require the use of a number of different R packages. Most of these can be sourced from CRAN, but a couple need to be downloaded from GitHub or BioConducter. We can use the following code to quickly load in the packages and install any packages not previously installed in the R console.

setwd("../data")

if (!require("pacman")) install.packages("pacman")
pacman::p_load("conover.test", "dendextend",  "ggdendro", "ggfortify", "multcompView","patchwork", "RColorBrewer", "Redmonder", "reshape2", "tidyverse", "vegan", "WGCNA")
pacman::p_load_gh("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")

Additional settings and color palettes

We’re also going to set the number of decimal places to display to 4 and create a few color palettes for our plots. colorRampPalette() creates a gradient between the colors specified (at least 2 colors are needed). By creating the object getPalette we can now call it and specify the number of colors we want to interpolate from the palettes we specify below. e.g. getPalette(10) will provide 10 hexadecimal color codes. Additionally we will make a color palette for our coral Family data where each of the four orders we are intesested in have a ramp of color in a different hue.

options("scipen" = 4)

colPal = divPal= c("#00B8AA", "#374649", "#FD625E", "#F2C811", "#788580")

getPalette = colorRampPalette(brewer.pal(8, "Dark2"))
getPalette2 = colorRampPalette(redmonder.pal(8, "qPBI"))

AL = colorRampPalette(c("#FFFFFF", getPalette2(4)[1]))(8)[3:8]
AN = colorRampPalette(c("#FFFFFF", getPalette2(4)[2]))(4)[2:4]
H = getPalette2(4)[3]
S = colorRampPalette(c("#FFFFFF", getPalette2(4)[4]))(11)[2:11]

coralPalette = c(AL,AN,H,S)

pie(rep(1,length(coralPalette)), col = coralPalette)

Import data

---

Coral count data

Import and prepare coral (Orders Scleractinia, Alcyonacea, Antipatharia, and Anthoathecata) and Scleractinian count data. These will be used for diversity metrics.

coralFamCounts = read.csv("coralFamilyCounts.csv")
coralFamCounts$Bank = factor(coralFamCounts$Bank, levels(coralFamCounts$Bank)[c(5,2,1,3,4)])
coralFamCountsBank = coralFamCounts %>% group_by(Bank) %>% summarise_each(sum)

sclerSppCounts = read.csv("scleractinianSppCounts.csv")
sclerSppCounts$Bank = factor(sclerSppCounts$Bank, levels(sclerSppCounts$Bank)[c(5,2,1,3,4)])
sclerSppCountsBank = sclerSppCounts %>% group_by(Bank) %>% summarise_each(sum)

Transformation of percentage and density by transect data

For our statistical analyses of transects by bank we will load in and transform the data from CPCe. We will standardize by transects and squareroot transform to deflate the influence of 0 counts which are commonly found in ecological abundance data.

mjrPerc = read.csv("nwgomMajorPerc.csv")
mjrPerc$Bank = factor(mjrPerc$Bank, levels(mjrPerc$Bank)[c(5,2,1,3,4)])
mjrPercT = mjrPerc
mjrPercT[3:ncol(mjrPerc)] = sqrt(mjrPerc[3:ncol(mjrPerc)] * 100)

mnrPerc = read.csv("nwgomMinorPerc.csv")
mnrPerc$Bank = factor(mnrPerc$Bank, levels(mnrPerc$Bank)[c(5,2,1,3,4)])
mnrPercT = mnrPerc
mnrPercT[3:ncol(mnrPerc)] = sqrt(mnrPerc[3:ncol(mnrPerc)] * 100)


mnrDens = read.csv("nwgomMinorDensity.csv")
mnrDens$Bank = factor(mnrDens$Bank, levels(mnrDens$Bank)[c(5,2,1,3,4)])
mnrDensT = mnrDens
mnrDensT[3:ncol(mnrDens)] = sqrt(decostand(mnrDens[3:ncol(mnrDens)],
                                                  "total") * 100)

Calculate percent cover and coral density by bank

We can group and calculate percent cover and density by bank and then rearrange them to make it easy to plot stacked barplots by bank. To do this we will use the melt() function in reshape2, which takes our “wide” data and melts it into “long” format.

mjrPercBank = mjrPerc[2:ncol(mjrPerc)] %>% group_by(Bank) %>% summarise_each(mean) %>%melt(id= "Bank", value.name = "PercentCover", variable.name = "MajorCategory")
mjrPercBank$MajorCategory = factor(mjrPercBank$MajorCategory, levels(mjrPercBank$MajorCategory)[c(1, 2, 9, 12, 4:7, 11, 3, 14, 10, 8, 13)])

coralPercBank = read.csv("coralFamilyPercentCoverBank.csv") %>% melt(id= "Bank", value.name = "PercentCover", variable.name = "Family")

coralPercBank$Bank = factor(coralPercBank$Bank, levels(coralPercBank$Bank)[c(5,2,1,3,4)])

coralPercBank$PercentCover = (coralPercBank$PercentCover) * 100

coralDensBank = read.csv("coralFamilyDensityBank.csv", check.names = FALSE, header = TRUE) %>% melt(id = "Bank", value.name = "Density", variable.name = "Family")

coralDensBank$Bank = factor(coralDensBank$Bank, levels(coralDensBank$Bank)[c(5,2,1,3,4)])

Now the data are in the correct format to build stacked barplots.

Percent cover and coral density by bank

---

First, we can examine the percent cover of major taxa/benthic category and the percent cover and density of corals by bank. We can make a series of barplots in R to see how each bank compares to the others.

Major category percent cover

colorCount1 = length(levels(mjrPercBank$MajorCategory))

mjrTypeA = ggplot(mjrPercBank, aes(x = Bank, y = PercentCover, fill = MajorCategory)) +
  geom_bar(stat = "identity", color = "black", size = 0.15) +
  scale_fill_manual(values = getPalette(colorCount1), name = "Major Category") +
  labs(y = "Percent Cover") +
  guides(fill = guide_legend(ncol = 2)) +
  theme_classic()

mjrType = mjrTypeA + theme(
  axis.title.x = element_blank(),
  axis.text.x = element_blank(),
  axis.title.y = element_text(color = "black", size = 8),
  axis.text.y = element_text(color = "black", size = 8),
  legend.position = "right",
  legend.justification = c(0,1),
  legend.title = element_text(color = "black", size = 8),
  legend.text = element_text(color = "black", size = 6),
  legend.key = element_blank(),
  legend.key.size = unit(0.5,"line"),
  legend.background = element_blank(),
  panel.background = element_blank(),
  panel.border = element_blank()
)

mjrType

Coral Family Percent cover plot

coralPercA = ggplot(coralPercBank, aes(x = Bank, y = PercentCover, fill = Family)) +
  geom_bar(stat = "identity", color = "black", size = 0.15) +
  scale_fill_manual(values = coralPalette, name = "Coral Family") +
  guides(fill = guide_legend(ncol = 2)) +
  labs(y = "Percent Cover") +
  theme_classic()

coralPerc = coralPercA + theme(
  axis.title.x = element_blank(),
  axis.text.x = element_blank(),
  axis.title.y = element_text(color = "black", size = 8),
  axis.text.y = element_text(color = "black", size = 8),
  legend.position = "none",
  legend.justification = c(0,1),
  legend.title = element_text(color = "black", size = 8),
  legend.text = element_text(color = "black", size = 6),
  legend.key = element_blank(),
  legend.key.size = unit(0.5,"line"),
  legend.background = element_blank(),
  panel.background = element_blank(),
  panel.border = element_blank()
)

coralPerc

Coral density bar plot

mjrDensA = ggplot(coralDensBank, aes(x = Bank, y = Density, fill = Family)) +
  geom_bar(stat = "identity", color = "black", size = 0.15) +
  # scale_fill_manual(values = densPalette, name = "Cnidarian Fam") +
  scale_fill_manual(values = coralPalette, name = "Coral Family") +
  guides(fill = guide_legend(ncol = 2)) +
  labs(y = expression(paste("Density (m"^"-2",")"))) +
  theme_classic()

mjrDens = mjrDensA + theme(
  axis.title.x = element_blank(),
  axis.text.x = element_text(color = "black", size = 8, angle = 45, vjust = 1,
                             hjust = 1),
  axis.title.y = element_text(color = "black", size = 8),
  axis.text.y = element_text(color = "black", size = 8),
  legend.position = "right",
  legend.justification = c(0,1),
  legend.title = element_text(color = "black", size = 8),
  legend.text = element_text(color = "black", size = 6),
  legend.key = element_blank(),
  legend.key.size = unit(0.5,"line"),
  legend.background = element_blank(),
  panel.background = element_blank(),
  panel.border = element_blank()
)

mjrDens

Group coral cover/density plots

Now that we have generated all 3 plots we can group them together with the package patchwork. This uses intuitive syntax to group plots, | separates plots next to each other and / separtes plots over one another. We can also gather similar legends into one legend and place it where we want. See the patchwork vignette for more details.

barPlots = (mjrType / coralPerc / mjrDens) +
  plot_annotation(tag_level = "A") &
  theme(plot.tag = element_text(size = 10, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_blank())

barPlots

ggsave("../figures/barPlots.eps", plot = barPlots, height = 12, width = 10,
       units = "cm", dpi = 300)

---

Alpha diversity metrics

---

Now we can examine the alpha diversity across different banks in the NW GOM. We can use vegan again to calculate different diversity metrics. We can then perform statistical tests to see where differences occur and create a figure to display the raw data that goes into these tests.

Coral Family alpha diversity metrics

First we will look at coral count data down to Family. We will build a new dataframe including richness, diversity, and evenness measures for all transects as well as banks as a whole.

coralFamRich = specnumber(coralFamCounts[,2:ncol(coralFamCounts)])
coralFamShannon = diversity(coralFamCounts[,2:ncol(coralFamCounts)],index = "shannon")
coralFamSimpson = diversity(coralFamCounts[,2:ncol(coralFamCounts)],index = "simpson")
coralFamRichBank = specnumber(coralFamCountsBank[,2:ncol(coralFamCountsBank)])
coralFamShannonBank = diversity(coralFamCountsBank[,2:ncol(coralFamCountsBank)],index = "shannon")
coralFamSimpsonBank = diversity(coralFamCountsBank[,2:ncol(coralFamCountsBank)],index = "simpson")

coralFamDiv = cbind(coralFamCounts[1],coralFamRich, coralFamShannon, coralFamSimpson)
colnames(coralFamDiv)[2:4] = c("rich", "shannon", "simpson")
coralFamDivBank = cbind(coralFamCountsBank[1], coralFamRichBank, coralFamShannonBank,
                        coralFamSimpsonBank)
colnames(coralFamDivBank)[2:4] = c("rich", "shannon", "simpson")

head(coralFamDiv)

##     Bank rich   shannon simpson
## 1 Bright    2 0.6931472   0.500
## 2 Bright    0 0.0000000   1.000
## 3 Bright    0 0.0000000   1.000
## 4 Bright    1 0.0000000   0.000
## 5 Bright    2 0.6931472   0.500
## 6 Bright    2 0.5623351   0.375

head(coralFamDivBank)

##       Bank rich  shannon   simpson
## 1 West FGB   17 1.759006 0.6972035
## 2 East FGB   18 2.270259 0.8612343
## 3   Bright    9 1.590340 0.7332906
## 4    Geyer    6 1.475558 0.7263522
## 5  McGrail   14 2.032716 0.8350812

Kruskal-Wallis test on coral Family diversity measurements

We can use a Kruskal-Wallis test on our non-normally distributed data to see how each metric differs across banks for the coral Family counts.

kWFamRich = kruskal.test(coralFamDiv$rich ~ coralFamDiv$Bank)

kWFamShannon = kruskal.test(coralFamDiv$shannon ~ coralFamDiv$Bank)

kWFamSimpson = kruskal.test(coralFamDiv$simpson ~ coralFamDiv$Bank)

kWFamRich$p.value

## [1] 0.00132383

kWFamShannon$p.value

## [1] 0.02010451

kWFamSimpson$p.value

## [1] 0.008953586

There are significant differences in coral Family Richness, Diversity, and Evenness across banks.

Conover-Iman test for pairwise comparisons

Following significant Kruskal-Wallis tests we can use Conover-Iman tests to perform pairwise comparisons among all banks. We will correct the p-values for multiple comparisons using FDR correction (Benjamini–Hochberg procedure).

richFamCIT = conover.test(x = coralFamDiv$rich, g = coralFamDiv$Bank, method = "bh", table = FALSE, list = TRUE, altp = TRUE)

##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 17.844, df = 4, p-value = 0
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -----------------------------------------
## Bright - East FGB   : -0.444526 (0.7300)
## Bright - Geyer      : -1.598914 (0.1849)
## East FGB - Geyer    : -1.681598 (0.1874)
## Bright - McGrail    : -2.439833 (0.0382)*
## East FGB - McGrail  : -2.986936 (0.0153)*
## Geyer - McGrail     : -0.810157 (0.5231)
## Bright - West FGB   :  0.302160 (0.7627)
## East FGB - West FGB :  1.278459 (0.2887)
## Geyer - West FGB    :  2.510978 (0.0419)*
## McGrail - West FGB  :  3.910724 (0.0011)*
## 
## alpha = 0.05
## Reject Ho if p <= alpha

shannonFamCIT = conover.test(x = coralFamDiv$shannon, g = coralFamDiv$Bank, method = "bh", table = FALSE, list = TRUE, altp = TRUE)

##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 11.6556, df = 4, p-value = 0.02
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -----------------------------------------
## Bright - East FGB   : -0.660928 (0.6365)
## Bright - Geyer      : -1.008450 (0.4487)
## East FGB - Geyer    : -0.626739 (0.5904)
## Bright - McGrail    : -2.137808 (0.1112)
## East FGB - McGrail  : -2.264480 (0.1214)
## Geyer - McGrail     : -1.173420 (0.4026)
## Bright - West FGB   :  0.170237 (0.8649)
## East FGB - West FGB :  1.420974 (0.3128)
## Geyer - West FGB    :  1.560563 (0.2993)
## McGrail - West FGB  :  3.304421 (0.0107)*
## 
## alpha = 0.05
## Reject Ho if p <= alpha

simpsonFamCIT = conover.test(x = coralFamDiv$simpson, g = coralFamDiv$Bank, method = "bh", table = FALSE, list = TRUE, altp = TRUE)

##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 13.5306, df = 4, p-value = 0.01
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -----------------------------------------
## Bright - East FGB   : -0.599600 (0.6865)
## Bright - Geyer      :  1.435765 (0.3043)
## East FGB - Geyer    :  2.648661 (0.0426)*
## Bright - McGrail    :  1.194780 (0.3886)
## East FGB - McGrail  :  2.524707 (0.0303)*
## Geyer - McGrail     : -0.374146 (0.7873)
## Bright - West FGB   : -0.707424 (0.6855)
## East FGB - West FGB : -0.193118 (0.8470)
## Geyer - West FGB    : -2.750835 (0.0632)*
## McGrail - West FGB  : -2.638367 (0.0293)*
## 
## alpha = 0.05
## Reject Ho if p <= alpha

Plot coral Family richness

First we can calculate group letters to represent the similarities between banks.

#set FALSE for n.s. comparisons and true for sig. comparisons
richFamDiff = ifelse(richFamCIT$altP.adjusted >= 0.05, FALSE, TRUE)

#name each comparison with the banks separated by only a dash, no space
names(richFamDiff) = richFamCIT$comparisons %>% str_replace(" - ", "-")

richFamSigLetters = multcompLetters(richFamDiff, compare="<", threshold=0.05,
                                     Letters = c(letters, LETTERS, "."), reversed = FALSE)
richFamLetters = as.data.frame(richFamSigLetters$Letters)
richFamLetters$Bank = row.names(richFamLetters)
colnames(richFamLetters)[1] = "id"

#set y for where letters will plot on y axis
richFamLetters$y = c(-0.95)

head(richFamLetters)

##          id     Bank     y
## Bright   ab   Bright -0.95
## East FGB ab East FGB -0.95
## Geyer    ac    Geyer -0.95
## McGrail   c  McGrail -0.95
## West FGB  b West FGB -0.95

Now we can add the letters to a violin plot including all transect points and a bar behind each bank to represent the bank total.

coralFamRichPlotA = ggplot(coralFamDiv, aes(x = Bank, y = rich, fill = Bank)) +
  geom_bar(data = coralFamDivBank, aes(x = Bank, y = rich, fill = Bank), stat = "identity", alpha = 0.6) +
  geom_violin(alpha = 0.5) +
  geom_point(shape = 21, color = "black",
             size = 1.5, position = position_jitterdodge(1)) +
  scale_alpha_manual(values = c(1, 1, 1, 1, 1)) +
  scale_fill_manual(values = divPal) +
  xlab("Bank") +
  ylab("Coral Family \nRichness") +
  annotate("label", label = paste("Kruskal-Wallis, p = ", round(kWFamRich$p.value, 4)),
           size = 3, x = 3, y = 19.8) +
  coord_cartesian( ylim = c(-0.79, 20)) +
  scale_y_continuous(breaks = seq(0, 20, 2))+
  geom_text(data = richFamLetters, aes(x = Bank, y = y, label = id),
            size = 4) +
  theme_classic()

coralFamRichPlot = coralFamRichPlotA +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.title.y = element_text(color = "black", size = 12),
        axis.text.y = element_text(color = "black", size = 10),
        legend.title = element_text(color = "black", size = 12),
        legend.text = element_text(color = "black", size = 10),
        panel.border = element_blank(),
        legend.position = "none"
  )

coralFamRichPlot

Plot coral Family diversity

# set FALSE for n.s. comparisons and true for sig. comparisons
shannonFamDiff = ifelse( shannonFamCIT$altP.adjusted >= 0.05, FALSE, TRUE)

# name each comparison with the banks separated by only a dash, no space
names(shannonFamDiff) = shannonFamCIT$comparisons %>% str_replace(" - ", "-")

shannonFamSigLetters = multcompLetters(shannonFamDiff, compare="<", threshold=0.05,
                                       Letters = c(letters, LETTERS, "."), reversed = FALSE)
shannonFamLetters = as.data.frame(shannonFamSigLetters$Letters)
shannonFamLetters$Bank = row.names(shannonFamLetters)
colnames(shannonFamLetters)[1] = "id"

#set y for where letters will plot on y axis
shannonFamLetters$y = c(-0.125)

head(shannonFamLetters)

##          id     Bank      y
## Bright   ab   Bright -0.125
## East FGB ab East FGB -0.125
## Geyer    ab    Geyer -0.125
## McGrail   a  McGrail -0.125
## West FGB  b West FGB -0.125

coralFamShannonPlotA = ggplot(coralFamDiv, aes(x = Bank, y = shannon)) +
  geom_bar(data = coralFamDivBank, aes(x = Bank, y = shannon, fill = Bank),
           stat = "identity", alpha = 0.6) +
  geom_violin(aes(fill = Bank), alpha = 0.5) +
  geom_point(aes(fill = Bank), shape = 21, color = "black",
             size = 1.5, position = position_jitterdodge(1)) +
  scale_fill_manual(values = divPal) +
  xlab("Bank") +
  ylab("Coral Family \nDiversity") +
  annotate("label", label = paste("Kruskal-Wallis, p = ", round(kWFamShannon$p.value, 4)),
           size = 3, x = 3, y = 2.48) +
  geom_text(data = shannonFamLetters, aes(x = Bank, y = y, label = id),
            size = 4) +
  scale_y_continuous(breaks = seq(0, 2.5, 0.5))+
  coord_cartesian(ylim = c(-0.105, 2.5)) +
  theme_classic()

coralFamShannonPlot = coralFamShannonPlotA +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.title.y = element_text(color = "black", size = 12),
        axis.text.y = element_text(color = "black", size = 10),
        legend.title = element_text(color = "black", size = 12),
        legend.text = element_text(color = "black", size = 10),
        panel.border = element_blank(),
        legend.position = "none"
  )

coralFamShannonPlot

Plot coral Family evenness

# set FALSE for n.s. comparisons and true for sig. comparisons
simpsonFamDiff = ifelse(simpsonFamCIT$altP.adjusted >= 0.05, FALSE, TRUE)

# name each comparison with the banks separated by only a dash, no space
names(simpsonFamDiff) = simpsonFamCIT$comparisons %>% str_replace(" - ", "-")

simpsonFamSigLetters = multcompLetters(simpsonFamDiff, compare="<", threshold=0.05,
                                       Letters = c(letters, LETTERS, "."), reversed = FALSE)
simpsonFamLetters = as.data.frame(simpsonFamSigLetters$Letters)
simpsonFamLetters$Bank = row.names(simpsonFamLetters)
colnames(simpsonFamLetters)[1] = "id"

#set y for where letters will plot on y axis
simpsonFamLetters$y = c(-.06)

head(simpsonFamLetters)

##           id     Bank     y
## Bright   abc   Bright -0.06
## East FGB   a East FGB -0.06
## Geyer     bc    Geyer -0.06
## McGrail    b  McGrail -0.06
## West FGB  ac West FGB -0.06

coralFamSimpsonPlotA = ggplot(coralFamDiv, aes(x = Bank, y = simpson, fill = Bank)) +
  geom_bar(data = coralFamDivBank, aes(x = Bank, y = simpson), stat = "identity",
           alpha = 0.6) +
  geom_violin( alpha = 0.5) +
  geom_point(shape = 21, color = "black",
             size = 1.5, position = position_jitterdodge(1)) +
  scale_fill_manual(values = divPal) +
  xlab("Bank") +
  ylab("Coral Family \nEvenness") +
  annotate("label", label = "Kruskal-Wallis, p < 0.0001 ",
           size = 3, x = 3, y = 1.19) +
  coord_cartesian( ylim = c(-.05, 1.2)) +
  geom_text(data = simpsonFamLetters, aes(x = Bank, y = y, label = id),
            size = 4) +
  theme_classic()

coralFamSimpsonPlot = coralFamSimpsonPlotA +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.title.y = element_text(color = "black", size = 12),
        axis.text.y = element_text(color = "black", size = 10),
        legend.title = element_text(color = "black", size = 12),
        legend.text = element_text(color = "black", size = 10),
        panel.border = element_blank(),
        legend.position = "none"
  )

coralFamSimpsonPlot

Scleractinina species alpha diversity metrics

sclerSppRich = specnumber(sclerSppCounts[,2:ncol(sclerSppCounts)])
sclerSppShannon = diversity(sclerSppCounts[,2:ncol(sclerSppCounts)],index = "shannon")
sclerSppSimpson = diversity(sclerSppCounts[,2:ncol(sclerSppCounts)],index = "simpson")

sclerSppRichBank = specnumber(sclerSppCountsBank[,2:ncol(sclerSppCountsBank)])
sclerSppShannonBank = diversity(sclerSppCountsBank[,2:ncol(sclerSppCountsBank)],index = "shannon")
sclerSppSimpsonBank = diversity(sclerSppCountsBank[,2:ncol(sclerSppCountsBank)],index = "simpson")

sclerSppDiv = cbind(sclerSppCounts[1],sclerSppRich, sclerSppShannon, sclerSppSimpson)
colnames(sclerSppDiv)[2:4] = c("rich", "shannon", "simpson")

sclerSppDivBank = cbind(sclerSppCountsBank[1], sclerSppRichBank, sclerSppShannonBank,
                        sclerSppSimpsonBank)
colnames(sclerSppDivBank)[2:4] = c("rich", "shannon", "simpson")

head(sclerSppDiv)

##     Bank rich shannon simpson
## 1 Bright    1       0       0
## 2 Bright    0       0       1
## 3 Bright    0       0       1
## 4 Bright    1       0       0
## 5 Bright    1       0       0
## 6 Bright    0       0       1

Kruskal-Wallis tests for Scleractinian spp. diversity metrics

kWSclerSppRich = kruskal.test(sclerSppDiv$rich ~ sclerSppDiv$Bank)

kWSclerSppShannon = kruskal.test(sclerSppDiv$shannon ~ sclerSppDiv$Bank)

kWSclerSppSimpson = kruskal.test(sclerSppDiv$simpson ~ sclerSppDiv$Bank)

kWSclerSppRich

## 
##  Kruskal-Wallis rank sum test
## 
## data:  sclerSppDiv$rich by sclerSppDiv$Bank
## Kruskal-Wallis chi-squared = 19.565, df = 4, p-value = 0.0006086

kWSclerSppShannon

## 
##  Kruskal-Wallis rank sum test
## 
## data:  sclerSppDiv$shannon by sclerSppDiv$Bank
## Kruskal-Wallis chi-squared = 2.4822, df = 4, p-value = 0.6478

kWSclerSppSimpson

## 
##  Kruskal-Wallis rank sum test
## 
## data:  sclerSppDiv$simpson by sclerSppDiv$Bank
## Kruskal-Wallis chi-squared = 27.513, df = 4, p-value = 0.00001566

Significant differences in Richness and Evenness of Scleractinian species among banks.

Conover-Iman test for pairwise comparisons

sclerSppRichCIT = conover.test(x = sclerSppDiv$rich, g = sclerSppDiv$Bank, method = "bh", table = FALSE, list = TRUE, altp = TRUE)

##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 19.5646, df = 4, p-value = 0
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -----------------------------------------
## Bright - East FGB   :  2.316291 (0.0425)*
## Bright - Geyer      :  1.828237 (0.1142)
## East FGB - Geyer    : -0.138445 (0.8900)
## Bright - McGrail    : -0.411932 (0.9724)
## East FGB - McGrail  : -3.611711 (0.0018)*
## Geyer - McGrail     : -2.614568 (0.0313)*
## Bright - West FGB   :  2.509386 (0.0316)*
## East FGB - West FGB :  0.361761 (0.7975)
## Geyer - West FGB    :  0.376388 (0.8836)
## McGrail - West FGB  :  3.838578 (0.0015)*
## 
## alpha = 0.05
## Reject Ho if p <= alpha

sclerSppSimpsonCIT = conover.test(x = sclerSppDiv$simpson, g = sclerSppDiv$Bank, method = "bh", table = FALSE, list = TRUE, altp = TRUE)

##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 27.5128, df = 4, p-value = 0
## 
## 
##                            Comparison of x by group                            
##                              (Benjamini-Hochberg)                              
## 
## List of pairwise comparisons: t statistic (adjusted p-value)
## -----------------------------------------
## Bright - East FGB   : -2.773601 (0.0118)*
## Bright - Geyer      : -1.911799 (0.0948)
## East FGB - Geyer    :  0.545434 (0.6510)
## Bright - McGrail    :  0.651177 (0.7364)
## East FGB - McGrail  :  4.555494 (0.0000)*
## Geyer - McGrail     :  2.977842 (0.0105)*
## Bright - West FGB   : -2.825232 (0.0126)*
## East FGB - West FGB : -0.124328 (0.9011)
## Geyer - West FGB    : -0.622405 (0.6677)
## McGrail - West FGB  : -4.592091 (0.0001)*
## 
## alpha = 0.05
## Reject Ho if p <= alpha

Scleractinian richness plot

#set FALSE for n.s. comparisons and true for sig. comparisons
sclerSppRichDiff = ifelse(sclerSppRichCIT$altP.adjusted >= 0.05, FALSE, TRUE)

#name each comparison with the banks separated by only a dash, no space
names(sclerSppRichDiff) = sclerSppRichCIT$comparisons %>% str_replace(" - ", "-")

sclerSppRichSigLetters = multcompLetters(sclerSppRichDiff, compare="<", threshold=0.05,
                                          Letters = c(letters, LETTERS, "."), reversed = FALSE)
sclerSppRichLetters = as.data.frame(sclerSppRichSigLetters$Letters)
sclerSppRichLetters$Bank = row.names(sclerSppRichLetters)
colnames(sclerSppRichLetters)[1] = "id"

#set y for where letters will plot on y axis
sclerSppRichLetters$y = c(-1.1)

head(sclerSppRichLetters)

##          id     Bank    y
## Bright   ab   Bright -1.1
## East FGB  c East FGB -1.1
## Geyer    ac    Geyer -1.1
## McGrail   b  McGrail -1.1
## West FGB  c West FGB -1.1

sclerSppRichPlotA = ggplot(sclerSppDiv, aes(x = Bank, y = rich, fill = Bank)) +
  geom_bar(data = sclerSppDivBank, aes(x = Bank, y = rich, fill = Bank),
           stat = "identity", alpha = 0.6) +
  geom_violin(alpha = 0.5) +
  geom_point(shape = 21, color = "black",
             size = 1.5, position = position_jitterdodge(1)) +
  scale_alpha_manual(values = c(1, 1, 1, 1, 1)) +
  scale_fill_manual(values = divPal) +
  xlab("Bank") +
  ylab("Scleractinian Species \nRichness") +
  annotate("label", label = paste("Kruskal-Wallis, p = ", round(kWSclerSppRich$p.value, 4)),
           size = 3, x = 3, y = 21.8) +
  coord_cartesian( ylim = c(-.9, 22)) +
  scale_y_continuous(breaks = seq(0, 22, 2))+
  geom_text(data = sclerSppRichLetters, aes(x = Bank, y = y, label = id),
            size = 4) +
  theme_classic()

sclerSppRichPlot = sclerSppRichPlotA +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(color = "black", size = 10, angle = 45, vjust = 1,
                                   hjust = 1),
        axis.title.y = element_text(color = "black", size = 12),
        axis.text.y = element_text(color = "black", size = 10),
        legend.title = element_text(color = "black", size = 12),
        legend.text = element_text(color = "black", size = 10),
        panel.border = element_blank(),
        legend.position = "none"
  )

sclerSppRichPlot

Scleractinian species diversity plot

sclerSppShannonPlotA = ggplot(sclerSppDiv, aes(x = Bank, y = shannon)) +
  geom_bar(data = sclerSppDivBank, aes(x = Bank, y = shannon, fill = Bank),
           stat = "identity", alpha = 0.6) +
  geom_violin(aes(fill = Bank), alpha = 0.5) +
  geom_point(aes(fill = Bank), shape = 21, color = "black",
             size = 1.5, position = position_jitterdodge(1)) +
  scale_alpha_manual(values = c(1, 1, 1, 1, 1)) +
  scale_fill_manual(values = divPal) +
  xlab("Bank") +
  ylab("Scleractinian Species \nDiversity") +
  annotate("label", label = paste("Kruskal-Wallis, p = ", round(kWSclerSppShannon$p.value, 4)),
           size = 3, x = 3, y = 2.48) +
  coord_cartesian( ylim = c(-0.1, 2.5)) +
  theme_classic()

sclerSppShannonPlot = sclerSppShannonPlotA +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(color = "black", size = 10, angle = 45,
                                   vjust = 1, hjust = 1),
        axis.title.y = element_text(color = "black", size = 12),
        axis.text.y = element_text(color = "black", size = 10),
        legend.title = element_text(color = "black", size = 12),
        legend.text = element_text(color = "black", size = 10),
        panel.border = element_blank(),
        legend.position = "none"
  )

sclerSppShannonPlot

Scleractinian evenness plot

#set FALSE for n.s. comparisons and true for sig. comparisons
sclerSppSimpsonDiff = ifelse(sclerSppSimpsonCIT$altP.adjusted >= 0.05, FALSE, TRUE)

#name each comparison with the banks separated by only a dash, no space
names(sclerSppSimpsonDiff) = sclerSppSimpsonCIT$comparisons %>% str_replace(" - ", "-")

simpsonSclerSppSigLetters = multcompLetters(sclerSppSimpsonDiff, compare="<", threshold=0.05,
                                            Letters = c(letters, LETTERS, "."), reversed = FALSE)
simpsonSclerSppLetters = as.data.frame(simpsonSclerSppSigLetters$Letters)
simpsonSclerSppLetters$Bank = row.names(simpsonSclerSppLetters)
colnames(simpsonSclerSppLetters)[1] = "id"

#set y for where letters will plot on y axis
simpsonSclerSppLetters$y = c(-.06)

head(simpsonSclerSppLetters)

##          id     Bank     y
## Bright   ab   Bright -0.06
## East FGB  c East FGB -0.06
## Geyer    ac    Geyer -0.06
## McGrail   b  McGrail -0.06
## West FGB  c West FGB -0.06

sclerSppSimpsonPlotA = ggplot(sclerSppDiv, aes(x = Bank, y = simpson, fill = Bank)) +
  geom_bar(data = sclerSppDivBank, aes(x = Bank, y = simpson), stat = "identity",
           alpha = 0.6) +
  geom_violin( alpha = 0.5) +
  geom_point(shape = 21, color = "black",
             size = 1.5, position = position_jitterdodge(1)) +
  scale_alpha_manual(values = c(1, 1, 1, 1, 1)) +
  scale_fill_manual(values = divPal) +
  xlab("Bank") +
  ylab("Scleractinian Species \nEvenness") +
  annotate("label", label = "Kruskal-Wallis, p < 0.0001 ",
           size = 3, x = 3, y = 1.19) +
  coord_cartesian( ylim = c(-.05, 1.2)) +
  geom_text(data = simpsonSclerSppLetters, aes(x = Bank, y = y, label = id),
            size = 4) +
  theme_classic()

sclerSppSimpsonPlot = sclerSppSimpsonPlotA +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(color = "black", size = 10, angle = 45, vjust = 1,
                                   hjust = 1),
        axis.title.y = element_text(color = "black", size = 12),
        axis.text.y = element_text(color = "black", size = 10),
        legend.title = element_text(color = "black", size = 12),
        legend.text = element_text(color = "black", size = 10),
        panel.border = element_blank(),
        legend.position = "none"
  )

sclerSppSimpsonPlot

Group all alpha diversity plots

diversityPlots = (coralFamRichPlot | coralFamShannonPlot | coralFamSimpsonPlot) / (sclerSppRichPlot | sclerSppShannonPlot | sclerSppSimpsonPlot) +
  plot_annotation(tag_levels = "A") &
  theme(plot.tag = element_text(size = 14, face = "bold"))

ggsave("../figures/sclerSppPlots.tiff", plot = diversityPlots, width = 20, height = 13,
       units = "cm", dpi = 300)

---

Hierarchical cluster analysis

---

We can visualize the dissimilarity and clustering of transects using hierarchical clustering analysis in R. We use the hclust() function with an agglomerative algorithm based on averages. We will generate dendrograms for the percent cover data (Major and minor taxa) and for coral density.

Major taxa/benthic substrate percent cover

mjrDend = mjrPercT[3:ncol(mjrPercT)] %>% vegdist %>% hclust(method = "average") %>% as.dendrogram

mjrDData = dendro_data(mjrDend)
mjrDendPoints = mjrDData$labels
mjrDendPoints$bank = mjrPercT[,2][order.dendrogram(mjrDend)]

mjrDendA = ggplot() +
  geom_segment(data = segment(mjrDData), aes(x = x, y = (1-y), xend = xend, yend = (1-yend)), size = 0.25) +
  geom_point(data = mjrDendPoints, aes(x = x, y = (1-y), fill = bank), size = 4, shape = 22, stroke = 0.25) +
  scale_y_reverse()+
  geom_hline(yintercept = 0.5, color = colPal[3], lty = 5) +
  scale_fill_manual(values = colPal, name = "Bank:") +
  guides(fill = guide_legend(override.aes = list(size = 4))) +
  coord_cartesian(xlim = c(10, 285)) +
  labs(y = "Similarity") +
  theme_classic()

mjrDend = mjrDendA + theme(
  axis.title.x = element_blank(),
  axis.text.x = element_blank(),
  axis.line.x = element_blank(),
  axis.ticks.x = element_blank(),
  axis.title.y = element_text(size = 14, color = "black", angle = 90),
  axis.text.y = element_text(size = 12, color = "black"),
  axis.line.y = element_line(),
  axis.ticks.y = element_line(),
  panel.grid = element_blank(),
  panel.border = element_blank(),
  panel.background = element_blank(),
  plot.background = element_blank(),
  legend.key = element_blank(),
  legend.title = element_text(size = 14),
  legend.text = element_text(size = 12))

mjrDend

Minor taxa percent cover

mnrDend = mnrPercT[3:ncol(mnrPercT)] %>% vegdist %>%
  hclust(method = "average") %>% as.dendrogram

mnrDData = dendro_data(mnrDend)
mnrDendPoints = mnrDData$labels
mnrDendPoints$bank = mnrPercT[,2][order.dendrogram(mnrDend)]

mnrDendA = ggplot() +
  geom_segment(data = segment(mnrDData), aes(x = x, y = (1-y), xend = xend, yend = (1-yend)), size = 0.25) +
  geom_point(data = mnrDendPoints, aes(x = x, y = (1-y), fill = bank), size = 4, shape = 22, stroke = 0.25) +
  scale_y_reverse() +
  geom_hline(yintercept = 0.4, color = colPal[3], lty = 5) +
  scale_fill_manual(values = colPal, name = "Bank:") +
  guides(fill = guide_legend(override.aes = list(size = 4))) +
  coord_cartesian(xlim = c(10, 285)) +
  labs(y = "Similarity") +
  theme_classic()

mnrDend = mnrDendA + theme(
  axis.title.x = element_blank(),
  axis.text.x = element_blank(),
  axis.line.x = element_blank(),
  axis.ticks.x = element_blank(),
  axis.title.y = element_text(size = 14, color = "black", angle = 90),
  axis.text.y = element_text(size = 12, color = "black"),
  axis.line.y = element_line(),
  axis.ticks.y = element_line(),
  panel.grid = element_blank(),
  panel.border = element_blank(),
  panel.background = element_blank(),
  plot.background = element_blank(),
  legend.key = element_blank(),
  legend.title = element_text(size = 14),
  legend.text = element_text(size = 12))

mnrDend

Coral density dendrogram

densDend = mnrDensT[3:ncol(mnrDensT)] %>% vegdist %>%
  hclust(method = "average") %>% as.dendrogram

densDData = dendro_data(densDend)
densDendPoints = densDData$labels
densDendPoints$bank = mnrDensT[,2][order.dendrogram(densDend)]

densDendA = ggplot() +
  geom_segment(data = segment(densDData), aes(x = x, y = (1 - y), xend = xend, yend = (1 - yend)), size = 0.25) +
  geom_point(data = densDendPoints, aes(x = x, y = (1 - y), fill = bank), size = 4, shape = 22, stroke = 0.25) +
  scale_y_reverse() +
  geom_hline(yintercept = 0.2, color = colPal[3], lty = 5) +
  scale_fill_manual(values = colPal, name = "Bank:") +
  guides(fill = guide_legend(override.aes = list(size = 4))) +
  coord_cartesian(xlim = c(10, 285)) +
  labs(y = "Similarity") +
  theme_classic()

densDend = densDendA + theme(
  axis.title.x = element_blank(),
  axis.text.x = element_blank(),
  axis.line.x = element_blank(),
  axis.ticks.x = element_blank(),
  axis.title.y = element_text(size = 14, color = "black", angle = 90),
  axis.text.y = element_text(size = 12, color = "black"),
  axis.line.y = element_line(),
  axis.ticks.y = element_line(),
  panel.grid = element_blank(),
  panel.border = element_blank(),
  panel.background = element_blank(),
  plot.background = element_blank(),
  legend.key = element_blank(),
  legend.title = element_text(size = 14),
  legend.text = element_text(size = 12))

densDend

Plot all dendrograms

percentDends = (mjrDend / mnrDend / densDend) +
  plot_annotation(tag_levels = "A") +
  plot_layout(guides = "collect") &
  theme(plot.tag = element_text(size = 16, face = "bold"),
        legend.position = "bottom", legend.background = element_rect(color = "black"),
  )

percentDends

ggsave("../figures/percentDends.eps", plot = percentDends, height = 15, width = 20, units = "cm", dpi = 300)

Non-metric multidimensional scaling

---

To look at community level differences among banks we can use non-metric multidimesional scaling (nMDS). With this ordination technique, points that are closer to onea another are more similar, and vice-versa.

Major Percent Cover nMDS

To run nMDS in R we can use the metaMDS() command in the package vegan. nMDS uses a (dis)similarity matrix, and we will use Bray-Curtis dissimilarity, which is the standard distance used in metaMDS(). Alternatively, if you wanted to use another dissimilarity measure you could specify another, i.e. metaMDS(dist, distance = "manahttan"). Since we already standardized and transformed the data we will set autotransform = FALSE. Since nMDS uses random restarts we can use set.seed() which allows randomized processes to be repeated with the same results if you use the same number for the standard input.

set.seed(694)
mjrNMDS = metaMDS(mjrPercT[3:ncol(mjrPercT)], trymax = 999, autotransform = FALSE, weakties = T)
mjrNMDS

mjrScores = as.data.frame(scores(mjrNMDS))
mjrScores$bank = factor(mjrPercT$Bank)
mjrScores$sample = row.names(mjrScores)
head(mjrScores)

##        NMDS1      NMDS2   bank sample
## 1 -0.7550395 -0.3394865 Bright      1
## 2 -0.6341406 -0.1007753 Bright      2
## 3 -0.5840860 -0.1380868 Bright      3
## 4 -0.6889981 -0.1849076 Bright      4
## 5 -0.6620902 -0.2192619 Bright      5
## 6 -0.8702873 -0.2859437 Bright      6

Build major percent nMDS biplot

After running the nMDS we use ggplot to create a biplot for visualizing the relationships between transects. Setting

mjrNMDSPlotA = ggplot() +
  geom_point(data = mjrScores, aes(x = NMDS1, y = NMDS2, fill = bank),
             color = "black", shape = 21, size = 2) +
  scale_fill_manual(values = colPal, name = "Bank:") +
  annotate("label", x = 0.32, y = 0.9, label = paste ("Stress = ", round(mjrNMDS$stress, 2), sep = ""), size = 4) +
  labs(x = "nMDS1", y = "nMDS2") +
  theme_bw()

mjrNMDSPlot = mjrNMDSPlotA +
  theme(axis.title.x = element_text(color = "black", size = 12),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.y = element_text(color = "black", size = 12),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "right",
        legend.title = element_text(color = "black", size = 12),
        legend.text = element_text(color = "black", size = 12),
        legend.key = element_blank(),
        legend.background = element_blank(),
        panel.border = element_rect(color = "black", size = 1.2),
        panel.background = element_rect(fill = "white"),
        panel.grid = element_blank(),
        plot.background = element_blank()
  )

mjrNMDSPlot

Minor percent cover nMDS

set.seed(694)
mnrNMDS = metaMDS(mnrPercT[3:ncol(mnrPercT)], try = 999, autotransform = F, weak.ties = T)
mnrNMDS

mnrScores = as.data.frame(scores(mnrNMDS))
mnrScores$bank = factor(mnrPercT$Bank)
mnrScores$sample = row.names(mnrScores)
head(mnrScores)

##        NMDS1      NMDS2   bank sample
## 1 -0.9138598 -0.5858688 Bright      1
## 2 -0.7282706 -0.5311320 Bright      2
## 3 -0.6650515 -0.4407828 Bright      3
## 4 -0.8045499 -0.2940012 Bright      4
## 5 -0.7685312 -0.5027324 Bright      5
## 6 -1.0598636 -0.1961360 Bright      6

Build minor percent cover nMDS biplot

mnrNMDSPlotA = ggplot() +
  geom_point(data = mnrScores, aes(x = NMDS1, y = NMDS2, fill = bank),
             color = "black", shape = 21, size = 2) +
  scale_fill_manual(values = colPal, name = "Bank:") +
  annotate("label", x = 0.3, y = 1, label =
             paste("Stress = ", round(mnrNMDS$stress, 2),
                   sep = ""), size = 4) +
  labs(x = "nMDS1", y = "nMDS2") +
  theme_bw()

mnrNMDSPlot = mnrNMDSPlotA +
  theme(axis.title.x = element_text(color = "black", size = 12),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.y = element_text(color = "black", size = 12),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "right",
        legend.title = element_text(color = "black", size = 12),
        legend.text = element_text(color = "black", size = 12),
        legend.key = element_blank(),
        legend.background = element_blank(),
        panel.border = element_rect(color = "black", size = 1.2),
        panel.background = element_rect(fill = "white"),
        panel.grid = element_blank(),
        plot.background = element_blank()
  )

mnrNMDSPlot

Coral dnsity nMDS

set.seed(694)
densNMDS = metaMDS(mnrDensT[3:ncol(mnrDensT)], trymax = 999, autotransform = F,
                   weak.ties = T)
densNMDS

densScores = as.data.frame(scores(densNMDS))
densScores$bank = factor(mnrDensT$Bank)
densScores$sample = row.names(densScores)
head(densScores)

##         NMDS1      NMDS2   bank sample
## 1 -0.38443677 -1.2796291 Bright      1
## 2 -0.36204672 -0.6250581 Bright      2
## 3 -0.81550294 -0.5951937 Bright      3
## 4 -0.09493783 -0.8846497 Bright      4
## 5 -0.15645331 -0.5426854 Bright      5
## 6  0.09616019 -1.4802012 Bright      6

Build coral density nMDS biplot

densNMDSPlotA = ggplot() +
  geom_point(data = densScores, aes(x = -(NMDS1), y = NMDS2, fill = bank),
             color = "black", shape = 21, size = 2) +
  scale_fill_manual(values = colPal, name = "Bank:") +
  annotate("label", x = 2.38, y = 1.68, label =
             paste("Stress = ", round(densNMDS$stress, 2),
                   sep = ""), size = 4) +
  labs(x = "nMDS1", y = "nMDS2") +

  theme_bw()

densNMDSPlot = densNMDSPlotA +
  theme(axis.title.x = element_text(color = "black", size = 12),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.y = element_text(color = "black", size = 12),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "right",
        legend.title = element_text(color = "black", size = 12),
        legend.text = element_text(color = "black", size = 12),
        legend.key = element_blank(),
        legend.background = element_blank(),
        panel.border = element_rect(color = "black", size = 1.2),
        panel.background = element_rect(fill = "white"),
        panel.grid = element_blank(),
        plot.background = element_blank()
  )

densNMDSPlot

Group all nMDS Plots

percNMDS = (mjrNMDSPlot | mnrNMDSPlot | densNMDSPlot) +
  plot_annotation(tag_levels = "A") +
  plot_layout(guides = "collect") &
  theme(plot.tag = element_text(size = 14, face = "bold"),
        legend.position = "bottom",
        legend.background = element_rect()
        )

percNMDS

ggsave("../figures/allNMDS.eps", plot = percNMDS, width = 20, height = 8, units = "cm", dpi = 300)

---

Permutational multivariate analysis of variance

---

nMDS and cluster analyses show some grouping of banks, with most differences appearing between E/WFGB and Geyer/Bright/McGrail. We can test for statiscally significant differences using PERMANOVA.

PERMANOVA

The adonis() function in vegan runs a PERMANOVA. We set Bank as a factor and use Bray-Curtis dissimilarity, as in our nMDS, and we will run 9,999 permutations. We will set the seed again, so that we get the same results here as in the publication.

set.seed(694)
mjrAdonis = adonis(mjrPercT[3:ncol(mjrPercT)] ~ Bank, data = mjrPercT, permutations = 9999, method = "bray")

set.seed(694)
mnrAdonis = adonis(mnrPercT[3:ncol(mnrPercT)] ~ Bank, data = mnrPercT, permutations = 9999, method = "bray")

set.seed(694)
mnrDensAdonis = adonis(mnrDensT[3:ncol(mnrDensT)] ~ Bank, data = mnrDensT, permutations = 9999, method = "bray")

mjrAdonis

## 
## Call:
## adonis(formula = mjrPercT[3:ncol(mjrPercT)] ~ Bank, data = mjrPercT,      permutations = 9999, method = "bray") 
## 
## Permutation: free
## Number of permutations: 9999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Bank        4     7.930 1.98259  21.264 0.22617 0.0001 ***
## Residuals 291    27.133 0.09324         0.77383           
## Total     295    35.063                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mnrAdonis

## 
## Call:
## adonis(formula = mnrPercT[3:ncol(mnrPercT)] ~ Bank, data = mnrPercT,      permutations = 9999, method = "bray") 
## 
## Permutation: free
## Number of permutations: 9999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Bank        4    16.627  4.1568  35.691 0.32913 0.0001 ***
## Residuals 291    33.892  0.1165         0.67087           
## Total     295    50.520                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mnrDensAdonis

## 
## Call:
## adonis(formula = mnrDensT[3:ncol(mnrDensT)] ~ Bank, data = mnrDensT,      permutations = 9999, method = "bray") 
## 
## Permutation: free
## Number of permutations: 9999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Bank        4    16.982  4.2455  13.756 0.18588 0.0001 ***
## Residuals 241    74.381  0.3086         0.81412           
## Total     245    91.363                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We can see that there are significant differences in both percent cover and coral density among banks.

Pairwise PERMANOVA

Since there were significant differences among banks, we can use pairwise PERMANOVA to perform bank by bank contrasts. Since we are doing so many pairwise comparisons we will correct the p-values using false discovery rate correction as we did in the Conover-Iman tests. We will run 9,999 permutations again.

set.seed(694)
mjrPWAdonis = pairwise.adonis(mjrPercT[3:ncol(mjrPercT)], factors = mjrPercT$Bank, sim.method = "bray", p.adjust.m = "BH", perm = 9999)

set.seed(694)
mnrPWAdonis = pairwise.adonis(mnrPercT[3:ncol(mnrPercT)], factors = mnrPercT$Bank, sim.method = "bray", p.adjust.m = "BH", perm = 9999)

set.seed(694)
mnrDensPWAdonis = pairwise.adonis(mnrDensT[3:ncol(mnrDensT)], factors = mnrDensT$Bank, sim.method = "bray", p.adjust.m = "BH", perm = 9999)

mjrPWAdonis

##                   pairs Df SumsOfSqs   F.Model          R2 p.value   p.adjusted sig
## 1    Bright vs East FGB  1 3.2708458 29.342200 0.188887499  0.0001 0.0001666667  **
## 2       Bright vs Geyer  1 0.1405261  4.898246 0.092597516  0.0114 0.0142500000   .
## 3     Bright vs McGrail  1 0.4460446  6.484060 0.097528043  0.0059 0.0084285714   *
## 4    Bright vs West FGB  1 3.8901519 48.723008 0.294002677  0.0001 0.0001666667  **
## 5     East FGB vs Geyer  1 2.5098831 23.601777 0.149755781  0.0001 0.0001666667  **
## 6   East FGB vs McGrail  1 1.8290633 15.708112 0.097138676  0.0001 0.0001666667  **
## 7  East FGB vs West FGB  1 0.1336446  1.221272 0.005980143  0.2693 0.2693000000    
## 8      Geyer vs McGrail  1 0.1321794  2.074891 0.029609618  0.1258 0.1397777778    
## 9     Geyer vs West FGB  1 3.2916232 43.102259 0.256404994  0.0001 0.0001666667  **
## 10  McGrail vs West FGB  1 2.6214511 29.206939 0.175726350  0.0001 0.0001666667  **

mnrPWAdonis

##                   pairs Df SumsOfSqs   F.Model          R2 p.value   p.adjusted sig
## 1    Bright vs East FGB  1 5.5523190 40.689923 0.244105475  0.0001 0.0001428571  **
## 2       Bright vs Geyer  1 0.3851128  5.885658 0.109224938  0.0001 0.0001428571  **
## 3     Bright vs McGrail  1 0.5248938  4.590287 0.071067753  0.0044 0.0055000000   *
## 4    Bright vs West FGB  1 5.8675890 57.167630 0.328233380  0.0001 0.0001428571  **
## 5     East FGB vs Geyer  1 6.6202886 52.826051 0.282755274  0.0001 0.0001428571  **
## 6   East FGB vs McGrail  1 5.4630420 38.882147 0.210307742  0.0001 0.0001428571  **
## 7  East FGB vs West FGB  1 0.1838229  1.454185 0.007112524  0.2119 0.2119000000    
## 8      Geyer vs McGrail  1 0.2803724  2.950887 0.041590555  0.0343 0.0381111111   .
## 9     Geyer vs West FGB  1 6.9788646 75.147362 0.375460166  0.0001 0.0001428571  **
## 10  McGrail vs West FGB  1 5.8160940 51.980888 0.275058967  0.0001 0.0001428571  **

mnrDensPWAdonis

##                   pairs Df SumsOfSqs   F.Model         R2 p.value   p.adjusted sig
## 1    Bright vs East FGB  1 5.5251413 18.067765 0.15174355  0.0001 0.0001666667  **
## 2       Bright vs Geyer  1 1.3038834  5.232629 0.10213470  0.0003 0.0004285714  **
## 3     Bright vs McGrail  1 0.6724264  2.372188 0.03995454  0.0343 0.0381111111   .
## 4    Bright vs West FGB  1 4.7120517 14.183160 0.13484250  0.0001 0.0001666667  **
## 5     East FGB vs Geyer  1 6.1356085 20.778646 0.15767840  0.0001 0.0001666667  **
## 6   East FGB vs McGrail  1 7.1254433 23.199076 0.15977427  0.0001 0.0001666667  **
## 7  East FGB vs West FGB  1 0.5635982  1.705656 0.01081544  0.1093 0.1093000000    
## 8      Geyer vs McGrail  1 0.7716780  2.864715 0.04100374  0.0131 0.0163750000   .
## 9     Geyer vs West FGB  1 5.9615775 18.743913 0.15653332  0.0001 0.0001666667  **
## 10  McGrail vs West FGB  1 6.4458338 19.607954 0.14898761  0.0001 0.0001666667  **